{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Mixed elements can be used to overcome locking when the material becomes incompressible. However, for the elements to be stable, they need to fulfil the LBB condition. We here show what happens with a linear / linear displacement pressure element (which does not fulfil the LBB condition). In the numerical example, we consider the Cook's Membrane problem with an applied traction on the right hand side."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using JuAFEM\n",
    "using BlockArrays"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function create_cook_grid(nx, ny)\n",
    "    dim = 2\n",
    "    corners = [Vec{dim}((0.0,   0.0)),\n",
    "               Vec{dim}((48.0, 44.0)),\n",
    "               Vec{dim}((48.0, 60.0)),\n",
    "               Vec{dim}((0.0,  44.0))]\n",
    "    grid = generate_grid(Triangle, (nx, ny), corners);\n",
    "    # Extract the left boundary\n",
    "    addfaceset!(grid, \"clamped\", x -> norm(x[1]) ≈ 0.0);\n",
    "    return grid\n",
    "end;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dim = 2\n",
    "# Interpolations\n",
    "\n",
    "ip_u = Lagrange{dim, RefTetrahedron, 1}()\n",
    "ip_p = Lagrange{dim, RefTetrahedron, 1}()\n",
    "\n",
    "# Quadrature rules\n",
    "qr           = QuadratureRule{dim  , RefTetrahedron}(3)\n",
    "qr_face      = QuadratureRule{dim-1, RefTetrahedron}(3)\n",
    "\n",
    "\n",
    "cellvalues_u = CellVectorValues(qr, ip_u);\n",
    "facevalues_u = FaceVectorValues(qr_face, ip_u);\n",
    "\n",
    "cellvalues_p = CellScalarValues(qr, ip_p);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# DofHandler\n",
    "function create_dofhandler(grid)\n",
    "    dh = DofHandler(grid)\n",
    "    push!(dh, :u, dim) # Add a displacement field\n",
    "    push!(dh, :p, 1)   # Add a pressure field\n",
    "    close!(dh)\n",
    "    return dh\n",
    "end;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Boundaryconditions\n",
    "function create_boundaryconditions(dh, grid)\n",
    "    dbc = ConstraintHandler(dh)\n",
    "    # Add a homogenoush boundary condition on the \"clamped\" edge\n",
    "    add!(dbc, Dirichlet(:u, getfaceset(grid, \"clamped\"), (x,t) -> zero(Vec{2}), [1,2]))\n",
    "    close!(dbc)\n",
    "    t = 0.0\n",
    "    update!(dbc, t)\n",
    "    return dbc\n",
    "end;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function symmetrize_lower!(K)\n",
    "    for i in 1:size(K,1)\n",
    "        for j in i+1:size(K,1)\n",
    "            K[i,j] = K[j,i]\n",
    "        end\n",
    "    end\n",
    "end;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "immutable LinearElasticity{T}\n",
    "    G::T\n",
    "    K::T\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function doassemble{dim}(cellvalues_u::CellVectorValues{dim}, cellvalues_p::CellScalarValues{dim},\n",
    "                         facevalues_u::FaceVectorValues{dim}, K::SparseMatrixCSC, grid::Grid,\n",
    "                         dh::DofHandler, mp::LinearElasticity)\n",
    "   \n",
    "    f = zeros(ndofs(dh))\n",
    "    assembler = start_assemble(K, f)\n",
    "    assembler2 = start_assemble(K, f)\n",
    "    nu = getnbasefunctions(cellvalues_u)\n",
    "    np = getnbasefunctions(cellvalues_p)\n",
    "    \n",
    "    global_dofs = zeros(Int, nu + np)\n",
    "\n",
    "    fe = PseudoBlockArray(zeros(nu + np), [nu, np]) # Local force vector\n",
    "    Ke = PseudoBlockArray(zeros(nu + np, nu + np), [nu, np], [nu, np]) # Local stiffness mastrix\n",
    "\n",
    "    t = Vec{2}((0.0, 1/16)) # Traction vector\n",
    "    ɛdev = [zero(SymmetricTensor{2, dim}) for i in 1:getnbasefunctions(cellvalues_u)]\n",
    "    for cell in CellIterator(dh)\n",
    "        fill!(Ke, 0)\n",
    "        fill!(fe, 0)\n",
    "        assemble_up!(Ke, fe, cell, cellvalues_u, cellvalues_p, facevalues_u, grid, mp, ɛdev, t)\n",
    "        celldofs!(global_dofs, cell)\n",
    "        assemble!(assembler, global_dofs, fe, Ke)\n",
    "    end\n",
    "    return K, f\n",
    "end;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function assemble_up!(Ke, fe, cell, cellvalues_u, cellvalues_p, facevalues_u, grid, mp, ɛdev, t)\n",
    "    \n",
    "    n_basefuncs_u = getnbasefunctions(cellvalues_u)\n",
    "    n_basefuncs_p = getnbasefunctions(cellvalues_p)\n",
    "    u▄, p▄ = 1, 2\n",
    "    reinit!(cellvalues_u, cell)\n",
    "    reinit!(cellvalues_p, cell)\n",
    "    \n",
    "    # We only assemble lower half triangle of the stiffness matrix and then symmetrize it.\n",
    "    @inbounds for q_point in 1:getnquadpoints(cellvalues_u)\n",
    "        for i in 1:n_basefuncs_u\n",
    "            ɛdev[i] = dev(symmetric(shape_gradient(cellvalues_u, q_point, i)))\n",
    "        end\n",
    "        dΩ = getdetJdV(cellvalues_u, q_point)\n",
    "        for i in 1:n_basefuncs_u\n",
    "            divδu = shape_divergence(cellvalues_u, q_point, i)\n",
    "            δu = shape_value(cellvalues_u, q_point, i)\n",
    "            for j in 1:i\n",
    "                Ke[BlockIndex((u▄, u▄), (i, j))] += 2 * mp.G * ɛdev[i] ⊡ ɛdev[j] * dΩ\n",
    "            end\n",
    "        end\n",
    "      \n",
    "        for i in 1:n_basefuncs_p\n",
    "            δp = shape_value(cellvalues_p, q_point, i)\n",
    "            for j in 1:n_basefuncs_u\n",
    "                divδu = shape_divergence(cellvalues_u, q_point, j)\n",
    "                Ke[BlockIndex((p▄, u▄), (i, j))] += -δp * divδu * dΩ\n",
    "            end\n",
    "            for j in 1:i\n",
    "                p = shape_value(cellvalues_p, q_point, j)\n",
    "                Ke[BlockIndex((p▄, p▄), (i, j))] += - 1/mp.K * δp * p * dΩ\n",
    "            end\n",
    "            \n",
    "        end\n",
    "    end\n",
    "    \n",
    "    symmetrize_lower!(Ke)\n",
    "\n",
    "    @inbounds for face in 1:nfaces(cell)\n",
    "        if onboundary(cell, face) && (JuAFEM.cellid(cell), face) ∈ getfaceset(grid, \"right\")\n",
    "            reinit!(facevalues_u, cell, face)\n",
    "            for q_point in 1:getnquadpoints(facevalues_u)\n",
    "                dΓ = getdetJdV(facevalues_u, q_point)\n",
    "                for i in 1:n_basefuncs_u\n",
    "                    δu = shape_value(facevalues_u, q_point, i)\n",
    "                    fe[i] += (δu ⋅ t) * dΓ\n",
    "                end\n",
    "            end\n",
    "        end\n",
    "    end\n",
    "end;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  1.188742 seconds (376.89 k allocations: 17.416 MiB, 0.59% gc time)\n",
      "  0.020581 seconds (10.04 k allocations: 298.641 KiB)\n"
     ]
    }
   ],
   "source": [
    "function solve(ν, doexport = true)\n",
    "    # Material\n",
    "    Emod = 1.\n",
    "    Gmod = Emod / 2(1 + ν)\n",
    "    Kmod = Emod * ν / ((1+ν) * (1-2ν))\n",
    "    mp = LinearElasticity(Gmod, Kmod)\n",
    "    # Grid, dofhandler, boundary condition\n",
    "    n = 50\n",
    "    grid = create_cook_grid(n, n)\n",
    "    dh = create_dofhandler(grid)\n",
    "    dbc = create_boundaryconditions(dh, grid)\n",
    "\n",
    "    # Assembly and solve\n",
    "    K = create_sparsity_pattern(dh);\n",
    "    @time K, f = doassemble(cellvalues_u, cellvalues_p, facevalues_u, K, grid, dh, mp);\n",
    "    apply!(K, f, dbc)\n",
    "    u = Symmetric(K) \\ f;\n",
    "\n",
    "    # Export\n",
    "    if doexport\n",
    "        vtkfile = vtk_grid(\"up_$ν\", dh)\n",
    "        vtk_point_data(vtkfile, dh, u)\n",
    "        vtk_save(vtkfile)\n",
    "    end\n",
    "    return u\n",
    "end\n",
    "\n",
    "for ν in [0.3, 0.4999999]\n",
    "    solve(ν)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compressible ν = 0.3\n",
    "![compressible.png](figures/mixed_up_compressible.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Incompressible ν  ≈ 0.5\n",
    "![incompressible.png](figures/mixed_up_incompressible.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.015089 seconds (10.04 k allocations: 298.641 KiB)\n",
      "Cook passed!\n"
     ]
    }
   ],
   "source": [
    "u = solve(0.3, false)\n",
    "Test.@test maximum(u) ≈ 26.13381519901358\n",
    "println(\"Cook passed!\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#= TODO: Mini element\n",
    "\n",
    "immutable MiniDisplacements{dim, shape, order} <: Interpolation{dim, shape, order} end\n",
    "JuAFEM.getnbasefunctions(::MiniDisplacements{2, RefTetrahedron, 1}) = 4\n",
    "\n",
    "function JuAFEM.value!(ip::MiniDisplacements{2, RefTetrahedron, 1}, N::AbstractVector, ξ::Vec{2})\n",
    "    @assert length(N) == 4\n",
    "    JuAFEM.value!(Lagrange{2, RefTetrahedron, 1}(), view(N, 1:3), ξ)\n",
    "    N[4] = N[1] * N[2] * N[3]\n",
    "    return N\n",
    "end\n",
    "\n",
    "function JuAFEM.derivative!{T}(ip::MiniDisplacements{2, RefTetrahedron, 1}, dN::AbstractVector{Vec{2, T}}, ξ::Vec{2, T})\n",
    "    @assert length(dN) == 4\n",
    "    ξx, ξy = ξ[1], ξ[2]\n",
    "    JuAFEM.derivative!(Lagrange{2, RefTetrahedron, 1}(), view(dN, 1:3), ξ)\n",
    "    dN[4] = Vec{2, T}((ξy * (1 - 2ξx - ξy),\n",
    "                       ξx * (1 - 2ξy - ξx)))\n",
    "    \n",
    "    return dN\n",
    "end\n",
    "\n",
    "cellvalues_u_mini = CellVectorValues(qr_mini, ip_u_mini);\n",
    "facevalues_u_mini = FaceVectorValues(qr_face_mini, ip_u_mini);\n",
    "ip_u_mini = MiniDisplacements{dim, RefTetrahedron, 1}()\n",
    "qr_mini      = QuadratureRule{dim  , RefTetrahedron}(3)\n",
    "qr_face_mini = QuadratureRule{dim-1, RefTetrahedron}(3)\n",
    "\n",
    "# Integrates along the right boundary\n",
    "function integrate_gamma(u, facevalues_u, grid, dh)\n",
    "    global_dofs = zeros(Int, ndofs_per_cell(dh))\n",
    "    u_integrated = zero(Vec{2})\n",
    "    for cell in CellIterator(dh)\n",
    "        celldofs!(global_dofs, cell)\n",
    "        up_nodes = u[global_dofs]\n",
    "        u_nodes = up_nodes[1:getnbasefunctions(facevalues_u)]\n",
    "        for face in 1:nfaces(cell)\n",
    "            if onboundary(cell, face) && (JuAFEM.cellid(cell), face) ∈ getfaceset(grid, \"right\")\n",
    "                reinit!(facevalues_u, cell, face)\n",
    "                for q_point in 1:getnquadpoints(facevalues_u)\n",
    "                    dΓ = getdetJdV(facevalues_u, q_point)\n",
    "                    u_cell = function_value(facevalues_u, q_point, u_nodes)\n",
    "                    u_integrated += u_cell * dΓ\n",
    "                end\n",
    "            end\n",
    "        end\n",
    "    end\n",
    "    return u_integrated\n",
    "end\n",
    "=#"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "name": "julia-0.6",
   "display_name": "Julia 0.6.0",
   "language": "julia"
  },
  "language_info": {
   "mimetype": "application/julia",
   "file_extension": ".jl",
   "version": "0.6.2",
   "name": "julia"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
